o 



43 



cd 



Path integral derivations of 
novel complex trajectory methods 



Jeremy SchifT 1 , Yair Goldfarb 2 and David J. Tannor 2 

o. 

Department of Mathematics, Bar-Ilan University, Ramat Gan, 52900 Israel 



2 
Department of Chemical Physics, The Weizmann Institute of Science, Rehovot, 76100 Israel 



July 30, 2008 



but taking into account the fi dependence of both the amplitude and the phase of the 
intial wave function, thus giving rise to the need for complex classical trajectories. 



£h ■ Abstract 

Path integral derivations are presented for two recently developed complex trajectory 

techniques for the propagation of wave packets, Complex WKB and BOMCA. Com- 
(N ■ 
£> ■ plex WKB is derived using a standard saddle point approximation of the path integral, 

ov 
in 

^1 , BOMCA is derived using a modification of the saddle point technique, in which the 

path integral is approximated by expanding around a near-classical path, chosen so 
that up to some predetermined order there is no need to add any correction terms to 
the leading order approximation. Both Complex WKB and BOMCA give the same 
leading order approximation; in Complex WKB higher accuracy is achieved by adding 
correction terms, while in BOMCA no additional terms are ever added — higher ac- 
curacy is achieved by changing the path along which the original approximation is 
computed. The path integral derivation of the methods explains the need to incorpo- 
rate contributions from more than one trajectory, as observed in previous numerical 
work. On the other hand, it emerges that the methods provide efficient schemes for 
computing the higher order terms in the asymptotic evaluation of path integrals. The 
understanding we develop of BOMCA suggests that there should exist near-classical 
trajectories that give exact quantum dynamical results when used in the computation 
of the path integral keeping just the leading order term. We also apply our path in- 
tegral techniques to give a compact derivation of the semiclassical approximation to 
the coherent state propagator. 



1 Introduction 

In a recent series of papers jTJ[2|,[3] we have considered two complex trajectory techniques for 
solving the time-dependent Schrodinger equation (TDSE). By a "trajectory technique" we 
mean that we solve the TDSE for the wave function by integrating a system of ODEs along 
certain trajectories in configuration space. By a "complex trajectory technique" we mean 
that the relevant trajectories evolve in complex configuration space -- i.e. we analytically 
continue the wave function and consider it as a function of complex space variables. (Note 
that the time variable remains real, so the trajectories are real curves in complex space.) 
The motivation for using complex trajectories comes from the substitution 

i> = exp (^£) , SeC, (1) 

in the TDSE 

h 2 

ihiPt = -— V 2 ^ + V{x)i> . (2) 

2m 

This yields the complex quantum Hamilton Jacobi equation (CQHJE) jH [5] 

S t +-±- (VS) 2 + \/(x) = ^-V 2 S . (3) 

2m 2m 

Taking h as small, the CQHJE can be considered as a perturbation of the classical Hamilton 
Jacobi equation (HJE) 

St+i L {vs)2 + v(x) = ° • (4) 

Since the classical HJE can be solved exactly by integration along trajectories in space 

defined by 

dx __ VS , . 

dt m 

it is natural to try a similar technique (at least as an approximation) for its perturbation, 

the CQHJE. Complex trajectories arise since S in equation (JTJ), and hence VS 1 in equation 

(J5J), are complex, leading to complex initial conditions for the evolution; furthermore the 

perturbation in ([3]) is complex. 

In our earlier papers we observed that a reasonable approximation to the wave function 
may require taking into account contributions from more than one trajectory reaching a 
particular point in space. We gave no theoretical justification for this, and one of the 
purposes of the current paper is to fill this hole. More generally the aim of the current paper 
is to strengthen the theoretical basis of the techniques of our previous papers by showing 
how they can be derived by saddle point evaluation of the wave function in the path integral 



representation. In this approach the need to (potentially) add the contributions of several 
trajectories emerges naturally. 

In section 2 we give a detailed presentation of the two techniques under study here. In 
the first method, Complex WKB [3], the trajectories are solutions of the classical equations 
of motion. In the second method, which we call BOMCA [HE], (BOhmian Mechanics with 
Complex Action), the trajectories are order % perturbations of solutions of the classical 
equations. Indeed, in BOMCA the trajectories depend on the order in h, while in complex 
WKB they remain the same; to increase the order in h in complex WKB we have to integrate 
a further system of ODEs along the trajectories. 

In section 3 we present a path integral derivation of Complex WKB. Since the trajectories 
involved in Complex WKB are classical, this involves a standard saddle point approximation 
of the path integral. The approach, however, is still nonstandard in that the initial value 
of the wave function is taken into account, leading to complex classical trajectories. In this 
regard, our approach differs from standard time-dependent WKB theory [Bj. Our approach 
is the appropriate one when the initial wave function has the form exp(iS imt (x)/h), as 
implied by ([I]) (though note an interesting recent paper of Maia et al. |7J). The equivalence 
of Complex WKB to the saddle point approximation of the path integral gives rise to a 
potentially very useful result. While it has long been recognized that certain factors involved 
in the (leading order) saddle point approximation of path integrals, specifically elements of 
the so-called stability matrix, can be calculated efficiently by integrating certain ODEs along 
classical trajectories [8], it turns out that the same is true (at least in our situation) for all 
higher order correction terms. Within the path integral formulation, the expressions for 
higher order correction terms involve complicated multiple integrals; the Complex WKB 
method reformulates these expressions as solutions of a system of first order ODEs, which 
is much easier to handle computationally. 

Before presenting the path integral derivation of BOMCA, in section 4 we present a 
slight modification to the standard method of asymptotic analysis of integrals with a large 
parameter. In section 5 we apply this modification to the path integral, and are led to 
BOMCA. The distinction between BOMCA and Complex WKB becomes extremely clear. 
Complex WKB and BOMCA give the same leading order approximation to the wave func- 
tion ^(X, T), determined as follows: First find complex classical trajectories x(t) satisfying 
appropriate initial and final conditions, specifically, solve the problem 

mx(t) + W(x(£)) = , x(0) = -— VlogVo(x(0)) , x(T) = X . (6) 

m 

Here ?/>o is the wave function at t — 0. Next, for each such trajectory, compute the matrix 



U satisfying 

mil + H(V)(*(t))U = , U(0) = I , 17(0) = --H (log Vo)(x(0)) . (7) 

m 

Here i7(l/) denotes the matrix of second derivatives of V and i7(log ipo) the matrix of second 
derivatives of log^o- Then the wave function is approximated by 

where S[x\ denotes the classical action associated with the path x(£), i.e. 

S[*] = [ (^x 2 - V(x)) (ft . (9) 

The sum in flH]) is over contributing trajectories (possibly not all trajectories), as we will 
explain later. The distinction between Complex WKB and BOMCA lies in the manner in 
which higher order corrections are made to (jHJ). In Complex WKB higher order corrections 
are made by multiplying the leading order contribution for each trajectory in (jSJ) by a 
suitable factor of the form 1 + 0(h). In BOMCA, the formula (jHJ) is never modified, but 
the paths x(t) and matrices U(t) are no longer required to be classical. More explicitly, the 
differential equations in (jBj) and (j7|) are replaced by equations of the form 

mx + W(x(i)) = 0(h) , (10) 

mU + H(V)(x(t))U = 0(h) . (11) 

BOMCA gives explicit expressions for the terms to introduce on the right hand side of these 
equations, but, as we shall explain, they are not unique choices. 

We call the quantity appearing on the right hand side of ([HD, with the classical choice 
of x and U, the classical wave function. Note that our use of the term "classical wave 
function" differs from previous uses, see for example Box 2.2 in [9]. We emphasize also that 
our classical wave function differs from the usual approximations made in time dependent 
WKB theory; the difference can be traced to different assumptions about the h dependence 
of the initial wave function, with our choice requiring the use of complex trajectories. 

As we have explained, BOMCA provides a prescription for making the formula (jHJ) more 
accurate, to any order in h, by changing the equations that x and U satisfy. We are led to 
conjecture that there may exist choices of x and U, satisfying (TTOl - tflTl) . such that formula 
dHJ) is exact. Unfortunately, at this stage we only know how to describe the right hand sides of 
equations (TTOT) and (TT7TT) perturbatively in h, and, as we have indicated above, there are many 
choices (one being associated with BOMCA). If there exist choices of x and U for which 



(151) is exact, the relevant trajectories x would be an interesting intermediate object between 
classical and quantum trajectories. The usual notion of quantum trajectories (in Bohmian 
mechanics) are the paths in (real) configuration space satisfying x = — Im (Vlog^(x, £)) 
(see the books jTU] and [H] for extensive discussion). One of the properties of these trajec- 
tories is that the velocity diverges at a node of the wave function, so near nodes quantum 
trajectories are qualitatively different from classical trajectories. In distinction to this, the 
non-classical trajectories that arise in BOMCA are always perturbations of classical tra- 
jectories. Certainly it is possible to express the wave function if) in the form (jSJ) only in 
certain regions of X, T space (like any other semiclassical formula, our formula suffers from 
problems related to caustics and Stokes' lines), but in these regions we conjecture that there 
exist non-classical trajectories which are perturbations of classical trajectories, which make 
the formula (jSJ) exact. 

After our derivation of BOMCA from the path integral, in section 6 we discuss the 
application of our ideas to the evaluation of other quantities in quantum mechanics. There 
is an extensive literature on the use of complex classical trajectories to compute the coherent 
state propagator, (see for example HH H21 H31 [131 HS1 [ISl HZl HBl HSl I2DI I2U ) , SLnd we show how 
some relevant formulae can be derived using our techniques. Section 7 contains concluding 
comments. An appendix provides the multidimensional derivation of the classical wave 
function JSJ); in the most of the main text we present the path integral derivations just in 
the one dimensional case. 

We conclude this introduction with a discussion of some relevant literature that has not 
yet been mentioned. The use of complex classical trajectories in semiclassical quantum me- 
chanics goes back to Stine and Marcus [22] and Miller and George [231 [231 125] , and numerous 
different applications have been subsequently presented (see for example [261 [27]). As far as 
we know, the first attempt to use complex classical trajectories to propagate wave packets 
is in the works of Huber, Heller and Littlejohn [281 [22] (the superposition of contributions 
from more than one trajectory also appears in this work). Our work, however, is closer 
to the rather different viewpoint of Boiron and Lombardi [30]. Other developments in this 
area include the extensive work of de Aguiar and collaborators [211 E21 [33], and very recent 
contributions of Chou and Wyatt [34J and Sanz and Miret-Artes [35]. There are a number 
of papers on time dependent WKB that we also found illuminating: [361 E3 EH [39] . Recent 
numerical work of Bender, Brody and Hook [3D], suggesting strong connections between 
complex classical dynamics and quantum dynamics, is very encouraging; such connections 
also provided the motivation for the detailed studies of complex classical dynamics of Kay 
and Shnerb [HI W2\ . Likewise, there are interesting connections between our current work 



and a series of papers by Poirier and collaborators [151 S31 SSI SSI El SB] • Poirier considers 
a decomposition of the wave function as a sum of (nodeless) terms; we suspect this decom- 
position is strongly linked, if not identical to, the decomposition implied by (jSJ) (at least in 
the time-dependent case [37]). Finally, we note that the hierarchy of ODEs in the BOMCA 
method can be viewed as a complex version of the DPM of Trahan, Hughes and Wyatt 
[39] which was originally developed for real trajectories (see the exposition in [9] and the 
comparison of BOMCA and DPM in [50]). 

2 The Complex WKB and BOMCA methods 

As explained in the introduction, the starting point for both the methods we consider in this 
paper is the substitution ip = e lS l n in the time-dependent Schrodinger equation to obtain 
the CQHJE ([3]). Here S(x,t) is complex. Both of our methods consist of approximating 
the CQHJE by a system of equations that can be solved by integrating along trajectories in 
complex configuration space. Both of our methods allow us to systematically improve the 
order of approximation in such a way that we might reasonably expect, in a suitable limit, 
to obtain exact results. 

The first method, complex WKB, proceeds by an expansion of S in powers of K. The 
relevant trajectories, irrespective of the order of approximation, are solutions of the complex 
classical equations of motion. Complex WKB is described in detail in subsection 2.1. The 
second method, BOMCA, involves a different approximation scheme that will be described 
in detail in section 2.2. The relevant trajectories depend on the order of approximation. 
The question arises as to whether these trajectories have a well-defined limit as the order 
of approximation is increased, and what this limit is. We will discuss this matter more in 
section 5. 

In subsection 2.3 we summarize some of the findings of our previous work on complex 
WKB and BOMCA that are relevant for understanding the rest of this paper. 

2.1 Complex WKB 

Writing 

oo 

s(x,t) = J2s n (*,t)h n (12) 

n=0 



and substituting in the CQHJE we obtain the following PDEs for the component functions 

S n (x,t): 

Sot + ^~ (VSo) 2 + V( X ) = , (13) 

2m 

V 9 1 n_1 ?' 

S n t + — -VS n = -— £ VS m • VS n _ m + — V 2 <S n _!, n = l,2,... . (14) 

m 2m £z x 2m 

Along with the TDSE we assume we are provided with an initial wave function ^o( x ) — 
"0(x, 0), and this provides us with an initial condition 5(x, 0) = —ih log ipo(x) for the 
CQHJE. How are we to choose appropriate initial conditions for the component functions 
S n l The archetypal form of the initial wave function we wish to consider is a nonnormalized 
Gaussian wave-packet, which in one spatial dimension takes the form 

/ (a + ia x )(x - x ) 2 ip (x - x Q )\ , . 

i/j {x) = exp I + I . (15) 

Here a ,a>i,x ,Xi are real constants, with a > 0, related to the expectations and variance 
of position and momentum via 

(x) = x , (p) = Po , (16) 

((x-x f) = ^, {{v-Vo?) = Ka \ +al) ■ (17) 

2ao 2clo 

For this choice of ipQ we have 

S(x, 0) = —ihlogipo(x) = i(ao + iai)(x — xo) 2 + po(x — xo) , (18) 

so it is natural to choose 

S o (x,0) = i(a + iai)(x -x ) 2 +p (x -x ) , (19) 

S n (x,0) = 0, n = l,2,.... (20) 

In greater generality, we assume that we can write the initial wave function as ^o( x ) = 
exp(iS init (x)/h) and take 5 (x,0) = S init (x) and 5 n (x, 0) = for n > 1. It is possible to 
generalize to the case that S*(x, 0) can be expanded in a Taylor series in h. 

In complex WKB we opt to solve the system of equations (TT3 l -( fT4l) by integrating along 

trajectories defined by 

dx = VSo 

dt m 

Writing v = ^^ and taking the gradient of f|T3|) gives 

dvi , „. 1 dV , , 



Thus along the trajectories we have 

^ = -ivy(x) . (23) 

dt m y ' y ' 

(Here j| denotes the Lagrangian derivative along the trajectories ^ = J^ + v ■ V.) We see 
the trajectories are simply classical trajectories. Note however that they are trajectories in 
complexified space. The initial condition for the TDSE gives the initial condition ^^(x) as 
a complex function of x, and thus 

v(0) = — V5 init (x(0)) (24) 

m 

is in general complex. 

To summarize to this stage: in complex WKB we choose to integrate along trajectories 
given by (|2J1 . and from ([TBI we deduce that these are actually classical trajectories, i.e. 

solutions of 

<^ x dv 1 „, ,, x /„„s 

-77 = v , — = W x , (25) 

at dt m 

with the complex initial condition (j2"4"l) . Also from f|T3|) we deduce that the evolution of So 
down these trajectories is given by 

db(j 1 2 



dt 2 



mv 2 - 7(x) . (26) 



To compute the evolution of Si, S2, . . . along the trajectories we need to use the equations 
( TT4"1) . We have written these equations with precisely the Lagrangian derivative of S n along 
the trajectories on the left hand side. Writing the first few equations out more explicitly we 
have 



dSi 



7T^V 2 So , (27) 



dt 2m 

f = -5S<™> , + ^" (28) 

^1 = _i VSl ■ VS 2 + -^-V 2 S 2 . (29) 

at m 2m 

We see that to find Si we need to follow the evolution of V 2 So along the trajectory, which 
may be found by calculating the second spatial derivatives of ( fT3l) . To find S2, we see 
from (I28p that we need to follow the evolution of first and second spatial derivatives of Si 
along the trajectories. These are obtained by taking two spatial derivatives of (J2T|) . but to 
integrate the resulting equations we also need third and fourth derivatives of So, obtained 
by further differentiation of ffTBl . 



D{n,d) 


n = 1 


n = 2 


n = 3 


general n 


d= 1 


4 


9 


16 


(n + 1) 2 


d = 2 


7 


22 


50 


l(n + l)(n + 2)(4n + 3) 


d = 3 


11 


46 


130 


|(n + l)(n + 2)(2n 2 + 6n + 3) 


general d 


\{d 2 + U + A) 


^ 4 + ...) 


Md" + ---) 


see (EU) 



Table 1: Total number of functions to be evolved along the trajectories as a function of 
dimensionality d and order n. The full expressions in the cases n = 2 and n = 3 with general 
dare i(d 4 +10d 3 +47d 2 +86d+72) and ^(d 6 + 21d 5 +205d 4 +1035d 3 +3034d 2 +4344d+2880) 
respectively. Note the numbers listed include a contribution of d for finding VSq = mv, 
which in practice is already determined when finding the trajectories. 

Proceeding in this manner we see that to obtain Sq, S\, . . . , S n we need to follow up to 
the 2n'th derivatives of So, up to the 2(n — l)'th derivatives of S\ etc. along the trajectory. 
The total number of derivatives of order up to i of a scalar function is 

d(d + l)...{d + i-l) (d + i)\ 



, d(d+l) d(d + l)(d + 2) 



+ 



(30) 



2 6 i\ d\i\ 

where d denotes the number of spatial dimensions. Thus the total number of functions we 
need to follow along the trajectories is given by 



D{n,d) 



^ jd + i) 
f^f dW 

l even 



(31) 



We are not aware of a closed form expression for this sum, but tabulate it in Table 1 for 
some low values of d and n. For fixed n and large d we note that D(n, d) ~ d 2n /(2n)\, i.e. 
the number of functions we need to follow increases polynomially with the dimension. 

At this juncture we write out in full the equations that we must integrate to obtain 
5*0, S±, 5 2 : To find the trajectories we solve Newton's equations ff25|) with the required initial 
condition f[2~4"|) . The gradient of So along the trajectories can be identified with mv, and we 
do not need to recompute it. Higher derivatives of S along the trajectories are determined 
by integrating the equations: 



dSi 



0,ik 



dt 

dSo t iki 
dt 

dSo : ikl m 

dt 



Vik 2^ So t ijSo % 



rn 



,jk j 



— Vikl 2^ {So y ijSojkl + So,kjSo,jli + SqjjSqj 



m 



Ik i 



(32) 
(33) 



Viklm 2-^i \PO,i3^0jklm + ^0,kj^>0,jlmi + ^O^j^OJmik + »- , 0,m//>->( 



m 



mj^OJikl) 



2^ \So,jklSo,jim + SojkmSojU + Soji m Sojikj 



m 



(34) 



Here V^ denotes 9 ^.q x etc, and S ^k denotes -£r§^r etc. Derivatives of S\ are determined 
by integrating the equations 

,,' — ~ — /_, So,jji /> So,ijS\j , (35) 

dt 2m ^-r 1 m ^ ' J ' J 

3 3 

— rf — = 7; — /^So,jjik /.So.ikjSij > (So.kjSnj ; H~ SoijSikj) ■ (36) 

dt 2m . m^ m^ 

3 3 3 

Finally So, Si, S 2 are obtained by integrating the equations (126|) . (1271) and (|28|) respectively. 
The initial conditions for all these equations are obtained from the initial condition for the 
TDSE via the function 5" mit (x). Explicitly, we have 

5 (0) = S init (x(0)) S ,i(0) = 5f lit (x(0)) S ,ij{0) = Slf (x(0)) ... 
Si(0) = Si,i(0) = ft«(o) = o ... 

S 2 (0) = %(0) = 5 2)ii (0) = 



Note that for an initial Gaussian wave packet (of the form ( fl5|) in one dimension) the 
function S mit is quadratic in the spatial variable, so only its first two derivatives will be 
nonzero. 

A few final notes before leaving our description of the complex WKB method: First, 
observe that when computing to order n in complex WKB we use the derivatives of the 
potential up to order 2n. We will see later how this emerges from the path integral approach. 
Second, observe that equation (J32l) is an example of a matrix Riccati equation [51j, and in 
particular, it has solutions that become infinite in finite time. These singularities are a 
manifestation of the phenomenon of caustics, which appear in almost every application of 
the WKB method. We note however that the singularities in the matrix Riccati equation are 
pole-type singularities, and it is possible, in a suitable sense, to integrate through them [52] . 
This is reminiscent of the fact that it is often possible to "regularize" caustics [531 IMl 1551 156] . 
We are currently investigating the singularity structure of the full system of equations (127|) - 
(p8|) . (l32|) - (!34|) . (!35l) - (l36j) [57] . Finally, we mention that although in this paper we work with 
the multidimensional Schrodinger equation in the form (j2j), assuming the mass matrix to 
be a multiple of the identity, there is no problem extending our formalism to work with a 
general positive definite mass matrix. 

2.2 BOMCA 

BOMCA is an alternate trajectory based approach for solving the CQHJE ([3]). Unlike 
complex WKB it does not involve an expansion in powers of K. Another distinction is 

10 



that in complex WKB the trajectories are classical paths, and in BOMCA they are not. 
Furthermore, the trajectories in BOMCA depend on the order of the approximation. 

In BOMCA we aim to integrate the CQHJE ([3]) by integrating along trajectories defined 

by 

dx VS 

—— = v , where v = . (38) 

dt m y J 

Differentiating the CQHJE we see that along these trajectories the velocity field v satisfies 

From the CQHJE we see that along such trajectories 

^ = ± v * _ V(x) _ ^v 2 S . (40) 

dt 2m K J 2m K J 

The problem integrating (|39|) and (140]) is that we have no information about the second 
and third derivatives of S that appear on the right hand sides. Borrowing an idea from 
complex WKB, we differentiate the CQHJE to find equations for the evolution of second 
and higher derivatives of S along the trajectories. At this stage we just write the equations 
for evolution of second, third and fourth derivatives: 

~df = ~ Vl] ~m^ SipSpi + 2m^ Sl]PP ' (41) 

"77 "~~ 'ijk / t \iJipiJpjk ~~1~ "JjpOpki "T O^pOpij J "T ~ / j ^ijkpp > \^^) 

(JjU lib rfl LillL ~. 

"77 -~ 'ijkl ~ / ^ iJjpiJpjkl i ^jp^pkli ~r iJkp^plij T" ^Ip^pijk) 

1 ih 

2_^ / [^'ijp^'pkl T" OikpOpjl + Oilpbpjk) + — 2_^i ^ijklpp ■ (,4o) 

lib „ Zttil „ 

Apparently things have not improved: on the right hand sides of these equations fifth and 
sixth derivatives of S appear. Now we can state the procedure of the BOMCA approxima- 
tion: in nth order BOMCA ignore all terms involving derivatives of S of order exceeding 
2n. Thus, in 1st order BOMCA the nonclassical term in (|39|) is taken to be zero and the 
trajectories are simply classical trajectories. The evolution f)40p for S, however, involves 
a nonclassical term with second derivatives; but these second derivatives are computed by 
integrating (J4"T1) down the trajectories, after ignoring the term with 4th order derivatives 
in (JH|). A comparison with the equations of complex WKB establishes that lowest order 
BOMCA is equivalent to lowest order complex WKB (i.e. complex WKB where only the 
terms So and Si are retained.) 



11 



Moving on to 2nd order BOMCA, a nonclassical term with third derivatives of S now 
remains in the equation for the trajectories (139]) . and fourth derivatives of S appear in (J4"T1) . 
The evolution of the third and fourth derivatives of S is given by fl42T )- fl43|) after ignorning the 
higher derivative terms. We observe that the resulting equations are precisely the same as 
equations (J3"3"l)-(134"1) that appeared in complex WKB. The trajectories, however, are different 
- thus 2nd order BOMCA is not equivalent to complex WKB of any order. The same is 
true for higher order BOMCA. Complex WKB and BOMCA however share the property 
that order n calculations involves derivatives of the potential V of order up to In. 

Note that ignoring the 5th and 6th derivative terms in (|42[) - (|43j) gives rise to order H 
errors in the 3rd and 4th derivatives of S. Through equations (|39|) and (14~TI) this gives rise 
to order h 2 errors in the trajectory x and the second derivative of S. At first glance it seems 
that the order h 2 error in x should give rise to an order h 2 error in S, as calculated from (1401) . 
But a careful calculation shows that the errors induced in S by both the error in x and the 
error in Sy are of order % , and thus we have achieved second order accuracy in S (and first 
order accuracy in S/H, which is what determines the wave function). A similar calculation 
shows that in nth order BOMCA, as described above, we achieve nth order accuracy in 
S. There is no evident benefit to truncating the BOMCA equations, say, by ignoring 4th 
derivatives but not 3rd. This point was not adequately appreciated in [H[2]. 

For clarity, we collect here the evolution equations for 2nd order BOMCA in the one 
dimensional case: 

S t = \mv 2 -V{x) + ^S" , (44) 



(45) 



dx 

Tt = *' 

^ = --V'(x) + ^S'", (46) 

dt m v ; 2m 2 v ; 

<7<?" 1 if) 

%. = -V"(x)--(S") 2 + —S"" , (47) 

d <?'" 3 

^_ = -V'"(x) - -S"S'" , (48) 

dt y ' m y ' 

d 1"" A 3 

^_ = -V'"\x) - -S"S"" - -{S'") 2 , (49) 

(Jju lib 'lb 

This system is a singular perturbation of the Newton's equations. The system can be 
somewhat simplified. Introducing a new variable f(t) defined (up to multiplication by a 
constant) by S" = y^-, we can solve the S'" and S"" evolution equations and find that the 
trajectories are determined by 

m^ = -V'{x{t)) + ^-S"\t), (50) 

12 



m % = -^mw® + h^jw { L - T f ^ { v ""^ u )) + ^ s '» 2 ) du ) ^ 

S'"(t) = j^(k- f Q V"\x{u))f{ufdu) . (52) 

Here K, L are constants of integration related to S'"(0), S""(0). 

We have not yet dealt with the question of initial conditions in BOMCA, but this is 
straightforward. Writing, as before, 

5 init (x) =-iftlogV>(x,0) , (53) 

we take 

S(0) = 5 init (x(0)) , (54) 

v,(0) = -Si nit (x(0)), (55) 

m 

Sij(0) = S\f (x(0)) etc. (56) 

The number of equations is easier to discuss in BOMCA than in complex WKB. In nth 

(d + 2n\ 
order BOMCA we retain derivatives of S up to order 2n, i.e. a total of functions. 

V d J 
We also need to integrate to find x; thus there are a total of 

( d )+* (57) 

functions. But this number cannot be directly compared with the number in complex WKB. 
In complex WKB first the trajectories are computed using Newton's equations. If the aim 
is to determine the wave function at X at time T we solve (1251) . with boundary conditions 
(j24p and x(T) = X. Once the trajectories are determined, the evolution of all the other 
functions along the trajectories is computed. In BOMCA it is necessary to solve for all the 
functions (admittedly a rather smaller number) in order to determine the trajectories; that 
is we look for a solution of the full BOMCA system satisfying initial conditions (|54p - (|56p 
and the final condition x(T) = X. At this stage we have not made a complete study of the 
relative efficiencies of the two approaches. 

2.3 The need for multiple trajectories 

We refer to our previous papers [1] , [2] , [3] for full details of the implementation of the methods 
and explicit numerical examples. In order to determine the wave function at position X and 
time T, it is necessary to find trajectories x(t) satisfying all the necessary initial conditions 
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and the condition x(T) = X. The missing initial data is simply the starting point of the 
trajectories, x(0). In every case we investigated we found there were multiple trajectories 
satisfying all the necessary conditions. That is, there are various possible choices of x(0), 
and we refer to the different possible choices as different "branches". In certain cases, 
the wave function associated with one branch gives an accurate result. In other cases, it 
is necessary to add the wave functions associated with more than one branch to get an 
accurate result. In still other cases, one branch gives an overwhelmingly large contribution 
and has to be discarded. For certain values of X and T there are transitions between the 
different behaviors, and in the neighborhoods of such transitions we could not get reasonable 
accuracy with our methods. 

The upshot of all this is that our derivation of the complex WKB and BOMCA equations 
starting from the CQHJE ([3]) is apparently not telling us the whole story. In the following 
sections we will present derivations of complex WKB and BOMCA starting from the path 
integral formulation of quantum mechanics. In this approach the existence of multiple 
branches, and the need to sometimes incorporate one, sometimes more, is easily explained. 
We also find a non-technical explanation of what the trajectories in BOMCA are. At the 
same time, we see that the equations we have presented in detail for complex WKB and 
BOMCA actually provide an efficient way to perform certain higher order perturbative 
calculations with path integrals. 

3 A path integral derivation of complex WKB 

The aim of this section is to show how complex WKB emerges from the path integral for- 
mulation of quantum mechanics. In Feynman's path integral formulation the wave function 
(we work for now in one space dimension) is written 



if)(X,T)= K(x ,X,T)ip (xo)dxo, (58) 

J — oo 

where ipo(x ) = ip(x , 0) is the initial wavefunction and 

K(x , X,T)=J Vx exp i^\ (59) 

is the propagator. The propagator is represented as a sum over all possible paths x(t), 
< t < T, satisfying the boundary conditions x(0) = xq and x(T) = X. S[x] denotes the 
classical action of the path x, given by 

S[x] = f -mx 2 (t) - V(x(t)) dt. (60) 

Jo 2 
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Inserting eg. (1591) into eg. (1581) . moving ipo{ x o) into the argument of the exponent and incor- 
porating the integration over x into / T>x yields an alternate version of the path integral 
formulation 

r (iS\x\ \ r fi(S[x] + S init (x(0)))\ 
4{X,T) = JVxexpl^+]ogiJ> (x{0))\=JVx&vl-±-^ % J \ , (61) 

where now / T>x represents now a sum over all possible paths satisfying the single boundary 
condition x(T) = X, and we have written, as before, S mit (x) = —ih log ijjo(x). 

The next step is to evaluate ip(X, T) using a saddle point approximation. To this end we 
consider the variation of the term in the exponential in the path integral, and in particular 
identify paths for which the first order variation vanishes. Replacing x by x + e in the term 
in the exponential we have 

S[x + e] + S' m]i (x(0) + e(0)) = F \m (x + if - V(x(t) + e(t)) dt + S in[t (x(0) + e(0)) , 

Jo 2 

= S[x] + S' mit (x(0)) + [ T (mxe - V'(x)e) dt + S init '(x(0))e(0) 

Jo 

+ [ (^™ 2 - \v'\xy) dt+~S^"(x(0))e(0r 

+ T A I 5 init(n) (a;(0))e(0)" - f V {n) (x)e n dt] . (62) 

After an integration by parts, and using the fact that e(T) = 0, as all paths have the same 
fixed end point, the linear terms in e become 

- /" (mx + V'{x)) edt + (S init '(x(0)) - mx(0J) e(0) . (63) 

Thus we deduce that in a saddle point approximation of l[61\) . the approximation will be a 
sum of contributions from classical paths satisfying the initial condition 



i(0) = -^'(*(0)) = --#^ . ( 64 ) 

These are exactly the complex classical paths that appear in complex WKB. 

Proceeding to look at the guadratic terms in (162]) . we want the variable over which we 
integrate in the path integral to be dimensionless, so we rescale e by writing 

e(t) = ^5(t) . (65) 

After this change the guadratic terms in (1^21) become 

hT (f (¥ 2 - L V "^^) * + is W, "W»)*(°) 2 ) • («) 
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We are now in a position to write the saddle-point approximation to (16"Tj) : 



<*P ( ! £ ^T (D T (5 MlW (x(0))i(0)» - l\^(x)S"(t)dtJj . 

Here the sum is over complex WKB paths, that is paths x(t) obeying the classical equations 
of motion and the initial condition (1M|) . However, as is usual in saddle-point approximations, 
more detailed calculations are necessary to decide which of these paths should be included 
in the sum. We will return to this point shortly. 

3.1 The lowest order approximation 

To compute the lowest order approximation we just need to evaluate the Gaussian integral 

/^exp («r (/; (i#(«) - i- n V(Atm>) * + ^3<o) 2 )) ■ (68) 

We recall that the integration here is over paths 5(t) obeying the single condition S(T) = 0. 
As usual, we compute this integral by dividing the interval [0,T] into N subintervals and 
discretizing. Appropriate (second order) discretization formulas for the various terms are 

T . N ( N-l N-l \ 

l S 2 (t)dt « y k 2 + 2 E $\ - 2 £ <to+ij , (69) 

[v"{x{t))8\t)dt » §(^"(x(0))5 2 +EV"(x(^))^ , (70) 

where <5j denotes S(iT/N). Using these, the discretized version of the path integral is 

Jd N Aexp(—AAA T ) (71) 

where A = ( 5 5\ 5 2 ... 5jv_i ), A denotes the tridiagonal N x N matrix 



/ 



.4 



q -1 

-1 2-^»(*(f)) -1 

-1 2-^V"(z(f 



(72) 
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and 

(At this point it is maybe worthwhile noting that for an initial Gaussian wavefunction, 
gmit (x(0)) has a positive imaginary part.) The measure d N A here includes a nontrivial 
iV-dependent normalization; it turns out this should be chosen so that 

(For the standard rules for Gaussian integrals see for example [60J; the correct choice of 
normalization is determined by checking that we get the correct result for a free particle.) 
The computation of the determinant det A proceeds as follows [61]: For n = 1,2, ... ,N, 
denote the determinant of the n x n matrix in the top left corner of A by D n . Then we have 

Z>, = , = l + rSm ""<r (0)) + O(iV-), (75) 



^-^-^HI))-!.^"?^!^^, m 



and for 3 < n < N 



^( , -^ v 'HTr]K- fl - (77) 

The recursion can be written in the equivalent form 

We need to determine det A = D^. The recursion and the initial conditions are such that as 
N — > oo, the D n will tend to samples of a function D(s), defined on the interval < s < T, 
obeying the differential equation D(s) = ——V"(x(s))D and initial conditions -D(O) = 1 and 
-D(O) = — S imt (a;(0)). The determinant we seek is simply det A = D(T). 

To summarize, we have arrived at the lowest order approximation for the contribution 
of the path x(t) in the sum ( 1671) : It is given by 

D(T) l h 



Here S[x] denotes the classical action associated with the path x(t), which is a solution 
of Newton's equations obeying the conditions x(T) = X and x(0) = — S imt (x(0)). The 
function 5* mit is determined by the initial wave function via 5' imt (a;) = — ih log ip(x, 0). The 
function D(s) is the solution of the initial value problem 

D(s) = -—V"(x(a))D , D(0) = 1 , D(0) = -^ init "(a;(0)) . (80) 

m m 
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We check that the above gives an exact result in the case of the free particle (V = 0) 
and initial Gaussian wave function 

ip{x,0) =exp I + I . (81) 

The initial wave function here has three parameters, Xq and po which are real and a which 
is complex, with positive real part. Classical paths take the form x(t) = A + Bt. The 
coefficients A, B should be determined by requiring 

X = A + BT, B = — (p + 2ia(A - x )) . (82) 

m 

The classical action along the path x(t) is then given by S[x] = \mB 2 T and D(T) = 
1 + ±S init "(x(0))T = 1 + 2^. Putting everything together we obtain 

«™ - w^-^+^^M^) (83) 



Til 



1 ,^l_ T ^'(x-« ) -m\*(x-«,->^ ■ ^ T 



III 



j _|_ 2iaT \ ^H + 2ioT\ ^ TU J % \ ' 111 J 2%m 



It is straightforward to check that this is the exact solution of the Schrodinger equation for 
the given initial condition. 

Having computed the lowest order approximation to the path integral in the one di- 
mensional case, we now state the generalization to the multidimemsional case, leaving the 
proof to an appendix. Once again the saddle point paths are exactly the trajectories that 
appeared in the complex WKB method, specifcially they are classical paths x(t) obeying 
the initial condition 

x(0) = -VS init (x(0)) , (84) 

m 

c.f. (1241) . as well as the final condition x(T) = X. The contribution from any such path to 
the wave function ip(X.,T) takes the form 

1 /z(S[x] + S init (x(0)))\ 

= exp — i '- . (85) 

D(T) \ h ) 

Here, as in the one-dimensional case, 5[x] denotes the classicial action associated with the 
path x(t). The factor D(T) is determined as follows: Denote by U(s) the d x d matrix 
solution of the initial value problem 

t)(s) = --H{V)(x{8))U , U(0) = I , 17(0) = -F(5 init )(x(0)) , (86) 

m m 
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where here H(V) and H(S imt ) denote the dx d matrices of second derivatives of V and S mit 
respectively. Then D{T) = det(U(T)). 

We now wish to compare the lowest order path integral results with the lowest order 
approximation in complex WKB in the previous section. The path integral results all appear 
in the paragraph above. For ease, we assemble here all the necessary equations from complex 
WKB. The trajectories are determined from 



dx dv 1 „,.„, s , „ s 



with boundary conditions 



v(0) = -VSJ* (x(0)) , x(T) = X . 

m 



The wave function is given by 



^(X,r)=exp[^^ + S 1 (T)) . (89) 



The evolution equations of the necessary quantities along the trajectories are 

-mv 2 - V(x) , (90) 



doo 1 2 



dt 2 

§ = lt Sa *> (91) 

^ = -V ik (Mt))--J2So,i j So, j k, (92) 

at m^ 

with initial conditions 

5 (0) = 5 init (x(0)) , ^(0) = Sf (x(0)) , 5i(0) = . (93) 

The correspondence is almost immediate. All that is necessary to do is to identify the matrix 
with entries Sqjj in complex WKB with the matrix product mUU~ l in the path integral 
approach. With this identification, the evolution equation (1921) coincides with the second 
order evolution equation (1861) for U. Also after this identification, the evolution equation for 
Si, ( I9TT) reads ^ = ^Ti^UU^ 1 ), with solution (taking into account the appropriate initial 



conditions) S\(t) = |logdet U(t), so e lSl = 1/JD(T), giving the prefactor in fl85|) . Finally, 
So in complex WKB is identified with 5[x] + 5 mit (x(0)) in the path integral approach. 

The path integral approach has added one significant piece of information over the direct 
complex WKB approach presented in the previous section. In the path integral approach we 
use the saddle point method for asymptotic evaluation of an integral. As is well known, when 

19 



there are multiple saddle points, it is sometimes necessary to take more than one into account 
to get an accurate approximation of the integral being studied. Deciding which saddle points 
contribute requires detailed analysis on a case-to-case basis. But at least we have found an 
explanation for the observations of our earlier work [2] , [3] that for certain values of X and T it 
is necessary to include the contibutions of multiple trajectories. (It is interesting to compare 
this explanation for the origin of multiple trajectories with that given by Miller [58], based 
on the implicit nature of the equations that generate dynamical canonical transformations. 
We suspect that Miller's explanation may correlate with the existence of multiple solutions 
of the classical HJE, a connection that would bring us full circle to an understanding of the 
need for multiple trajectories in Complex WKB and BOMCA.) In future work [59] we hope 
to study the possible criteria for demarking different regions in X, T space in which different 
(numbers of) trajectories contribute. This is strongly interconnected with the existence of 
caustics. Caustics are points X,T at which the determinant D{T) vanishes (on at least one 
trajectory, in fact such points are associated with coallescing trajectories). It is possible to 
study the dynamics of such points, and from this to deduce certain information about the 
dynamics of the regions in which different numbers of trajectories contribute. Unfortunately, 
however, at the moment deciding on which trajectories to include in a calculation is more 
of an art than a science. 

3.2 The first order correction 

We now consider the higher order terms in f|67|) . We restrict ourselves in this section to the 
1-dimensional case. The series in the exponential in the third line of (loTj) is an expansion 
in half-integer and integer powers of the dimensionless parameter ^2 where L denotes a 
typical length scale of the functions V(x) and S mit (x). We are assuming this parameter 
is small. All terms with half-integer powers multiply odd powers of 5 and thus do not 
contribute to the value of the integral. The lowest order correction terms arise when we 
replace the exponential by 

hT ( ginit"/, Q r T 



(x(0))<5(0) 3 - f V"'(x)5 3 (t)dt\ . (94) 



72m 3 y Jo J 

After discretizing, in this approximation the integral ( 17TT) is replaced by an expression of 
the form 

N-l N-1N-1 



J d N A exp (^ AAA T ) 1 + £ 8} A + £ £ 5*$% , (95) 



HN 

iAA 1 j [ 1 + 

8=0 8=0 j=0 
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where Ai and Bj,-, which do not depend on the components of S, are 



^S initW, (x(0)) + 0U) 2 = 



a = { um Zs,sr;^: KNj . : , w 



ihT 3 

' 24m 2 N 



V""(x(§)) z>0 



B %0 



2 



?&^(o)r (f )) + o O) i > o, j = o 
^^ w (x(o))y- (x (f )) + o (£) i = o,; > o 

&^(§))v w (*(§)) z,i>o 



(97) 



HT° 



72m 3 N 2 



-A-J -> D(U)D(t s ) / -^, (100) 



The Gaussian integrals in the above expression are standard. Taking into account our 
normalization of the measure d n A we have 

/d"Aexp(fA^)^| = -_|=(3A^M-.M- + 2(A S 1 ) 3 ) (99) 

(c.f. [60J). A calculation similar to the calculation of det A in the previous subsection shows 
that as iV — > oo 

^A". 1 -► DHADCtA [ T -$- 

N l ° Kl) K J) J m , x ( tl ,t 3 ) D(u) 

where D(s) is the solution of fl86l) . (This calculation uses the fact that a second, linearly 
independent, solution of the differential equation in fl86|) is given by D(s) J S ^ .) In this 
manner we can write down the first order approximation to the path integral. For simplicity 
we restrict ourselves here to the case that 5* mit and S imt vanish, as otherwise the relevant 
formulae are lengthy. Combining the formulae above we find that in this case the first order 
approximation is found by multiplying the leading order approximation by 

+ ^3 [ [ v '" M*i)) v '" M«) D(hfD(kf 

du \ ( r T du \ ( f T du 




(ti,t 2 ) D{u) 2 J \Jt\ D{u) 2 J \Jt 2 D(u) 2 I Wmax(ti,t 2 ) D(u 



We now need to compare this with a similar term arising in the complex WKB method. 
The first order approximation in complex WKB is obtained by multiplying the leading order 
approximation by e %hS2 ^ T \ To find S2 it is necessary to integrate 5 new differential equations 
along the trajectories (in addition to those that have to be solved to find the lowest order 
approximation): the equations for Sq, Sq", S[, S" and S 2 , equations (J551) . (j53|) . (j55|) . (J551) and 
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respectively. All of these are linear equations, and an explicit formula can be written 
for the answer. To simplify matters we assume the initial conditions for all the 5 relevant 
quantities are zero, which is consistent with the assumption made in writing fllOll) . In its 
most obvious form (without making any attempts to simplify the integrals that appear) the 
solution takes the form 

+ *b f im (/»' ^w or ""wo)^)*) *•) 2 * 

' / ^ V'"(x(w))D 3 (w)dw )dv) du)dt . 



o D 2 {v) \Jo 

It is a straightforward but tedious matter to check that the factor fllOip is equal to 1 + ihSi 
(the first order approximation to e 2 ). 

Thus we see explicitly the equivalence of the first order approximation to the path integral 
and results from the complex WKB method retaining terms up to order S^. For consistency, 
this equivalence must continue to higher orders. Note that if we keep terms up to order h n in 
the path integral the resulting formulae will involve derivatives of V (and S imt ) up to order 
2n + 2, and the same is true if we retain terms up to S n+ \ in complex WKB. We note that 
in practice, complex WKB is far easier to implement for higher order corrections. Although 
the number of differential equations that need to be integrated along the trajectories grows 
rapidly with the order, as described in the previous section, it remains relatively easy to write 
down the necessary differential equations, and integrating the relevant first order system 
along the trajectories is easily handled using standard computer packages. Direct application 
of path integral methods involves the calculation of iterated integrals, as in (j!02p or fllOll) . 
which is a less standard procedure. The coefficients of the different iterated integrals (the 
number of which grows rapidly as order increases) also involve tricky combinatoric factors. 



4 A modification of standard asymptotic analysis 

In the previous section we have explained the connection of the complex WKB method as 
described in section 2 and the standard asymptotic evaluation of the path integral. We 
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would like to also understand BOMCA from this viewpoint. But there is a clear problem 
- whereas the trajectories in complex WKB are classical paths, corresponding to minima 
of the classical action, the paths in BOMCA are nonclassical. How can nonclassical paths 
possibly arise in the context of an asymptotic evaluation of the path integral? In this section 
we describe a modification of standard asymptotic analysis for Laplace-type integrals. In 
the next section we will apply what we have learn here to path integrals. 

The usual approach to asymptotic evaluation of integrals such as Jf^ g(x)e~ x ^ x ^dx, 
where A is a large positive parameter, proceeds as follows: The integral is dominated by 
contributions from regions close to the minima of f(x). Sufficiently near a minimum xq the 
function f(x) is approximated by a quadratic Taylor polynomial f(x ) + ^f"(xo)(x — xq) 2 . 
So we rewrite the integral in the form 

^e-^^+^o)^-^) 2 )^ , (103) 



oo 



where g(x) = g(x) exp (— A (f(x) — f(x ) — ^f"(x )(x — x ) 2 )), and evaluate the contribu- 
tion from the region near x by expanding g(x) in a Taylor series in x — Xq and evaluating 
the resulting integrals exactly. This gives a series in negative powers of A. 

The modification to this procedure that we want to consider is as follows: The Taylor 
polynomial approximation to f(x) at its minimum is only one of many ways to approximate 
f(x) in the appropriate region by a quadratic function. Suppose we choose another quadratic 
approximant. How does this change the resulting asymptotic expansion? 

For definiteness, we consider a specific example, asymptotic approximation of the facto- 
rial function for large n using the integral representation 

/•oo 

n\ = / e nlogx - x dx . (104) 

Jo 

The function in the exponent has a minimum at x — n, and the usual asymptotic formula 
for n! is obtained by approximating this function by the quadratic nlogn — n — ^(x — n) 2 
and rewriting the integral 

n\ ~ n n e~ n f°° e- {x ~ n ^' 2n ~g(x) dx , (105) 

J — oo 

where 

g(x) = exp ( nlogx — x — (nlogn — n) -\ (x — n) 2 

V 2n 

= l + ^(x-nf-^-nr + ^(x-nf + ^(x-nT + ...(106) 
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Integrating gives the standard asymptotic series for n\ 

, 1 / 1 1 13Q \ 

n\~V2^n n+ *e- n (l + + —- + ... . (107) 

V Yin 288n 2 51840™ 3 / y ' 

Note that to get the correct coefficient of n~ r it is necessary to keep certain terms of order 
up to 6r in the Taylor series ( 11061) . 

Suppose now that instead of using the above quadratic approximant for the exponent we 
use the more general approximant nlogiV — N — j§(x — N) 2 . Here S and N are currently 
undetermined, but for definiteness we assume that N = n + 0(1) and S = n + 0(1). The 
integral now takes the form 



rv. 



N n e -N f°° exp f / _ N y\ g( x )dx (108) 



«p ( (£ - <* - "> + \ (I - w) <* - ">' + t^r 1 ^ - N r 



where 

g(x) = exp ( n log x - x - (n log iV - iV) H — -(x - N) 2 ) (109) 

1/1 n\, „,„ y^ ( .,„_! n 

r=3 

Making the substitution x — N = \f§y this becomes 

, Am -n ^ r -vV2 ((n-N)VS N 2 -nS 2 ™. ^^nS^ 2 \ 
n\~N n e N VS]_J y/2 exp\^ ^ y + ^ y 2 + E(~l) r *— j/'J dy . 

(110) 
Note that in the second exponent here the coefficients of y and y 3 behave as n -1 / 2 , the 
coefficients of y 2 and y A behaves as n~ l , and in general for r > 3 the coefficient of y r 
behaves as n 1_T-//2 . We compute the integral by expanding the second exponential term 
in a power series in y and computing the resulting integrals exactly. The leading order 
approximation is y2ii~SN n e~ N . The first order correction arises from replacing the second 
exponential by 

, (N 2 -nS , nS 2 ,\ l(f(n-N)VS\ nS 3 ' 2 3 \ 2 ,,„, 

and integrating, to obtain 

, n N ( 1 / QN 6 + 10n 2 S 3 - QN 4 nS + QN 4 Sn 2 - 12N 5 Sn- 

n\ ~ y/2^SN n e~ N 1 + — — 

^ 12AT6 ^ 6Ar 6 5 _ 9^2^ + i2n 2 S 2 N 2 - 12nS 2 N 3 

' (H2) 

It can be verified directly that provided N = n + 0(1) and S = n + 0(1) the correction 

term is ^-- + 0(??r 2 ) for large n. Expanding the second exponential in (11 101) to suitable 
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higher order gives us higher correction terms, apparently depending on N and S as well as 
n, but in fact independent of the choice of N, S to the desired order. 

Essentially what we have shown above is that Stirling's formula for n! can be made to 
depend on two variables N and S while retaining all its properties. The obvious question 
that needs to be asked at this stage is whether N and S can be chosen usefully. A full 
investigation of this would take us off on a tangent to the main topic of this paper, so we 
limit ourselves here to the simple observation which will allow us to give a path integral 
derivation of BOMCA: It is possible to choose N and S as functions of n in such a way 
that all correction terms to the leading order approximation n\ ~ \ / 2itSN n e~ N vanish. 
Furthermore, at least in the case of the factorial function that we are looking at now, this is 
not simply a perturbative result; that is, we can find analytic functions N and S of n, with 



the correct asymptotic behavior for large \n\ and such that T(n + 1) = y / 2irSN n e~ N at least 
in some region of the complex plane including the positive real axis. We will see in the next 
section how BOMCA is related to an analogous result for path integrals. Presumably there 
should be some way to select N and S "well" on the basis of properties of the integrand of 
(1104[) . but we do not attempt to study this here. 

5 A derivation of BOMCA from the path integral 

The path integral is 

^(X,T) = |DxexpQ(5[x]+5 init (x(0)))) , (113) 

where as before the integration is over all paths with x(T) = X. Applying the idea presented 
in the previous section means approximating S[x] + S mit (x(0)) with a quadratic, which we 
will take of the form 

S[X] + S^(X(0)) + £ l -m{x{t) - X(t)f - l - (V"(X(t)) + q(t)) (x(t) - X{t)fdt 

+ ±(S^"(X(0)) + q(0))(x(0) - A(0)) 2 . (114) 

Here X(t) is the path around which we are expanding, still to be fully determined, but 
assumed to be an order % perturbation of a classical path. Likewise the function q(t) (which 
plays the role of S in the previous section) is currently undetermined, assumed of order h. 
Using this quadratic as our leading order approximation in the path integral gives 



MX,T) = ex P Q(5[X] + 5 init (X(0))))|De 
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exp I - U me{tf - (V"(X(t)) + q(t)) e{tfdt + (S^"(X(0)) + <z(0))e(0)< 
exp Nr (/ mX(t)e(t) - V(X(t))e(t) + ^M*) 2 " E TT^W * 



r=3 

oo cinit( r )/ 



^ init '(X(0)) e (0) - g(0) e (0) 2 + £ S " Ut r f m e(0) r j j • (115) 

Here we have written e(t) = x(t) — X (t) . We can simplify the second exponential in the path 

■ -+(r) 

integral in two ways. First, purely for ease of presentation we will assume that S = 

for r > 3, i.e. that the initial wave function is a Gaussian wave packet. There is no difficulty 
to restore the extra terms, but the calculations become extremely lengthy. Second, we make 
choices on the initial values of the currently unknown functions X(t) and q(t) to eliminate 
other terms in the second exponential as follows: First, we assume q(0) = 0. Second, 
integrating the term / mX(t)e(t) gives a boundary contribution — mX(0)e(0) and we can 
cancel this by requiring mX(0) = S mit (X(0)). Implementing all these simplifications gives 
us 

i/,(X,T) = ex P Q-(S[X] + S init (X(0)))) j^e (116) 

exp (± ( jT me(t) 2 ~ (V"(X(t)) + q(t)) e{tfdt + S init "(X(0)) e (0) 2 ) ) 
exp (± (j*(-mX(t) - V'(X(t)Mt) + l -q{t)e{tf - £ ^f ^ W eft)) . 

Finally, we move to dimensionless quantities by substituting e(t) = J^6(t), giving 

rl>(X,T) = expQ(5[X] + 5 init (X(0))))|^ 

exp ( % I- (fm 2 - - (V"(X(t)) + q(t)) 8{tfdt + -S^"(X(0))5(0) 2 )) 
\ 2 \Jo m m J J 




T ImT /_••_,., 1 „„„,.^ »,., T .. ,, 



h 



X(t) + -V'(X(t))) S(t) + —q(t)5(tf 
m J 2m 



. E !!W2%,|,. (iit) 



., m i r\ 



Assuming that both X(t) + —V'(X(t)) and q(t) are of order h, we see that the coefficients of 
6(t) and 5(t) 3 in the second exponential are of order h~ ' , the coefficients of 5(t) 2 and <5(£) 4 
are of order h , and in general the coefficient of 5{t) r is of order Ti r ' _ for r > 3. (This is in 
direct analogy to the caclulations for the factorial function in section 4.) The leading order 
approximation to the path integral is obtained by simply discarding the second exponential 
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term. The remaining Gaussian integral is identical to one we have already computed (but 
with V"(X(t)) replaced by V"(X(t)) +q(t)), and we obtain the leading order approximation 



exp {^(S[X} + S mit (X(0))) 
MX,T)= KnK [ \ ' (118) 



where f(s) is the solution of 



1 .„„.„. ,, , „ „, , .,„, . S ,„, 1 



" , 



f(s) = (V"(X(s)) + q(s))f(s) , /(0) = 1 , /(0) = -S^"(X(0)) . (119) 

m m 

To obtain the first order correction to this, we need to replace the second exponential in the 
path integral by 



2m \Jo h 12m 



n, nT ( f r (*(«) + ±v<*M» + irw dt 



h 6m 2 

There are 5 terms here that we need to consider, as oppposed to 2 in the derivation of the 
first correction term in Complex WKB. In addition to the integration formulae (I98l - (l99l we 
need the formulae 

/^exp(fA^)^ = jij;-^*, (121) 

I d " A »* (t aa&t ) ®> - -irTsn* 1 *? ■ (122) 

Computing all the necessary integrals gives the following result for the first order correction: 
The leading order approximation should be multiplied by 

1 - L£«ww([m) dt (m 

im r T r T „ T/ . „ w , „. , „. . / r T du 



- S / / NMNMmfih) I f , ^tt\ dhdt 2 

2n JO JO \Jmax(ti,t 2 ) J{U) Z J 

+ iff ^)v-W*,))/(*.)/(«' (/l (tife) ^(jf ^) ** 

"'' '' 'V" (*(«,)) v""Wt 2 ))/(*i) 3 /(*2) 3 



24m 3 jo jo 

3 (/ T *)(r *)(r *) +2 (r *Vi«, 

\Jmax(ti,t 2 ) j{uYJ \Jti JKUy ) \Jt 2 /(«) / \Jmax(ii,t 2 ) J (li) 
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Here we have written N(t) = X(t) + — V'(X(t)). We have left the integrals here in the 
form they arise using the relevant rules for Gaussian integrals. To manipulate the integrals, 
though, it is more convenient to write them in terms of integrals in which all the variables 
are all ordered. Doing this gives: 

1 - si;f *>£* i^f^jw* (124) 

7777 pT fT rT i 

dh / dt 2 / dh N(h)f(h)N(t 2 )f(t 2 " 



2 
J V"3V 

1 f T f T f T f T 



n Jo Jti Jt 2 j(t 3 ) 



+ - f dh f dt 2 I dh I dh N{h)f{h)V'"{h)f{hf-——— 
m Jo Jt! Jt 2 Jt 3 J{h> J{h) 

1 r T 1 r T^ T r T~ I 1 "1 

+ — ( dh I dh I dh I dh V"Xt 1 )f(t 1 ) 3 - F7 -r=N(t 3 )f(t3 



2 



2mA, L h *Jt 2 °A 3 * v LU v ±7 f(h) 2 w/JW/ /(t 4 ) 2 

'' dh I dt 2 I dt ?> I dh V"\h)f{hfN{h)f{h) 



mJo Ai Jta A 3 ' /O3) 2 f(h) 

+ -^ [ T dh [ T dt 2 [ T dh V""(h)f(h) 4 ' ' 



2 



2 



4m 2 A, V tl 'A 2 ° V1/JV ^ f(h) 2 f(h) 

+ -^ f T dh f T dt 2 f T dt 3 f T dh f T dt 5 v"'{h)f{hf^—,v'"{h)f{hT ' 



2m 3 A> V tl 'A 2 % "A* J ^ /JVi/ /(i 2 ) 2 vo/JVO/ f(h) 2 f(t 5 y 

^\iT) fT fT rT rT rT 1 I "1 

+ — - / dh dt 2 / dt 3 / dh / dh V"'{h)f{hfV"\h)f{hf 



2m* Jo L J tl z Jt 2 d A 3 *J U ° V1 " V1 ' v /JV 7 f(h) 2 f(h) 2 f(h) 2 

Our intention now is to choose the functions iV(i) and g(t) (both assumed to be of order 
ft) in such a way that there is no first order correction, i.e. so that the sum of the integrals 
in the above expression vanishes. From the above we see immediately that for any choice 
of N(t) it is possible to choose q(t) such that the first order correction terms vanish. One 
choice that suggests itself for N(t) is simply to take N(t) = 0. Then the correct choice of 
q(t) is 

q(t) = ^(5i' d " /( " )V "" (x( "» (125) 

1 rt ru rv 1 

+- du dv dwV , "{X{u))f{uf—-V"\X{w))f{wf 
m Jo Jo Jo J\ v ) 

+A f du r dv r dw i v'"(x(v))f(v) 3 v"\x(w))f(w) 3 ) . 

m Jo Jo Jo j{u) J 

The solution that is of main interest for us, however, is 

*<*> = -2^i V " ( - Y( " ))/( '" )3<fa < 126 > 

*> = ^(ii'* /( " )4 "'""(- Y( "» < 127 > 

Q ft ru rv I \ 

+- du dv dw-—V'"{X{v))f{vYV'"{X{w))f{wY) 
m Jo Jo Jo j{u) z J 
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We summarize what we have shown up to this point. For either of the choices of N(t) 
and q(t) given above (or for any other choice of N(t) and the appropriate matching choice 
of q(t)) we have demonstrated that the leading order approximation to the path integral, 
(IllSp . requires no first order correction. Here the path X and the function / are chosen to 
satisfy 

X(t) + —V'(X(t)) = N(t) , mX(0) = S init '(X(0)) , X(T) = X , (128) 

m 

mf(s) + V"(X(s))f(s) = -q(s)f(s) , /(0) = 1 , /(0) = ^ init "(X(0)) . (129) 

It is straightforward to check that the choice ( 1126l) -( TT27l) describes BOMCA (compare 
equations (HJ^ . (TT26I) . (TT27h . fTT28|) . f[T29T) with (jIlft.(150l),(lgTl).(152]): the constants K,L should 
be chosen to be zero, and recall that in the discussion of BOMCA we wrote S" = *j-i)- 
We have arrived at the understanding of BOMCA set out in the introduction -- that it 
corresponds to an evaluation of the path integral around a near-classical path, chosen in 
such a way that the classical wave function remains accurate to any desired order in h, 
with the path x and U being modified appropriately. We have found that in fact there are 
other ways to change x and U in such a way as to "correct" the classical wave function. 
In particular, we can continue to use classical paths, but replace the usual Jacobi equation 
with f)129p . where q is given by (J125J) . (The existence of this option extends to higher 
dimensions.) In practice the direct derivation of BOMCA, as given in section 2, is clearly 
preferable over the path integral for determining higher order corrections. The path integral 
approach, however, is necessary to understand the need to add contributions from different 
(near-) classical trajectories can be justified. 

Both the direct approach to BOMCA and the path integral approach only allow us to 
construct the relevant near-classical trajectories order-by-order in h. The question arises 
as to whether it is possible to find pairs x and U for which the classical wave function 
is exact. As we have already explained in the introduction, the relevant paths would be 
an intermediate object between classical paths and the quantum trajectories of Bohmian 
mechanics --on the one hand the new paths would be order h perturbations of classical 
paths, but on the other hand they would enable permit the derivation of exact quantum 
dynamical results, at least in regions of configuration space where they exist. At this stage 
the existence of such paths remains just a conjecture. 
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6 The coherent state propagator 

This section is a slight digression from the main point of this paper, but provides another 
illustration of our path integral methods, as well as the need for complex classical trajectories 
in "semiclassical" calculations. The so-called coherent state propagator has been studied 
extensively by many authors [HI H21 [131 HH HS1 [161 HZl HH1 UHl I2QI [2TJ . By the coherent 
state propagator we mean the overlap between the wave function ip(x, T) evolving from an 
(initial) coherent state with another (final) coherent state. In fact our methods allow us to 
write down a rather more general object; we write down the leading order approximation 
for the overlap between the wave function ?/>( x ? T) evolving from an initial state of form 
exp(iSi(x)/h) with a final state of the form exp(iSf(x)/h). Using the Feynmann path 
integral representation of the (standard) propagator, the overlap takes the form 

V = |_" &c/V/(x/) J™ dxiiptixi) J V* exp ( ^M J (130) 

where the path integral is over all paths satisfying x(0) = x$, x(T) = x/. We can absorb 

the integrations over the initial and final position into the path integral to write this simply 

as 

(g[x] + ft(x(0)) - S}(x(T))) 



V =-- I Dx exp I -* '- | (131) 

where now the path integration is over all paths x(t), < t < T, with no specified boundary 
conditions. 

To compute the semiclassical approximation to V we replace x in the exponent in the 
above expression by x + e and expand to second order in e. We choose x so that the 
linear term vanishes. Taking the action to be given by J |mx 2 — V(x) dt we find that the 
appropriate paths must satisfy 

mx + W(x) = 0, (132) 

mx(0) = VSi(x(0)) , (133) 

mk(T) = V5J(x(T)) . (134) 

The leading order approximation is thus a sum over such paths of the form 

?«^^(x(T))e X p(i5[x]/ft)^(x(0))JPe exp(Q[e]) (135) 

where 

Q[e] = ^([me 2 -e{t) T H{V)^{t))e{t)dt 

+e{Q) T H{S i ) (x(0))e(0) - e(T) T H(S})(x(T))e(T^j . (136) 
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Following the method of the calculation in appendix A, the factor JT>e exp(Q[s]) can be 
replaced (modulo some normalization factor) by the large N limit of 1/vdet V, where 

f I + Gq -I ... \ 

-I 21 + Gi -/ 

-/ 2I + G 2 -I 

1"= -I 2I + G 3 -I . (137) 

-/ 27 + GW-i -/ 

-I I + G N J 



Here 



Gq 
G r 

Gn 



T 



■fT(5i)(x(0)) 



r 2 



m jy--v-'/v--v-// 2mN 2 
T 2 ( (rT\ 

T 



ff(V)(x(0)) 



iV-1 



tf(S»(x(T)) 



r- 



;H(V)(x(T)) 



(138) 
(139) 
(140) 



The method for evaluating the determinant used in appendix A (with a slight addition to 
understand the nontrivial normalization) yields the final result 



V 



'2nih\ 
, m ) 



d/2 



E 



^(x(r))exp(zff[x]/ft)^(x(0)) 



\ 



det 



U X {T) + it/ 2 (T)ir(5i)(x(0)) - ^(fiy)(x(r))^(T)- 
-4/f(5t)(x(T))C/ 2 (T)i?(5 i )(x(0)) 



Here £/i and C/ 2 are two solutions of the equation 

ml7(t) = -H(V)U(t) , 
satisfying initial conditions 

C/i(0) = I f C/ 2 (0) = 

LTi(O) = ' \ J7 2 (0) = / 

From equation (11431) it follows that the entries of U2 have dimensions of time, whereas those 
of U\ are dimensionless. 

Restricting to the case of Gaussian initial and final states, taken in the form 

m(x-x j) T ^i(x-x i) ipoi ■ (x-Xoj)\ 



(141) 
(142) 

(143) 



^/(x) 



cxp 
exp 



2h h J 

m(x-x /) T fi/(x-x /) ip /-(x-x /) 



2/i 



ft 



(144) 
(145) 
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the above formula reduces to 

V « ( 2lTih ) ^ V ^( x ( T )) ex P(^[ x ]/^)^( x (°)) (146) 

'" ' ^ /det (U^T) + tU 2 {T)tti + XVfU^T) + Q. f U 2 {T)Q. % ) 



Here x j, poi are real parameters giving the expectation values of the position and momentum 
in the intial state, x /, po/ are real parameters giving the expectation of the position and 
momentum in the final state, and fl$, Qf are symmetric, complex matrices (with eigenvalues 
with positive real part). The relevant paths in the case of Gaussian initial states are those 
satisfying the boundary conditions 

rax(O) = poi + zmfij(x(0) - Xoj) , (147) 

m±(T) = p 0/ - imO}(x(T) - x 0/ ) , (148) 

c.f. HU H21 H31 H31 HSl HSl [TZl HSl UHl EQl I3T] . Note that the formula dHHD is not dimensionless 
as for simplicity we have been working with nonnormalized Gaussian states. 

We give one further simplifcation, just as an illustration of the use of this formula: In 
the scalar case for a free particle (Ui(t) = 1, Uzify = t) the formula gives the exact result 



V 



exp 



(149) 



{Pi ~ Pf) 2 + m 2 QiQ}{xi - x f f + iT{pfQ} + p 2 f Q { ) + 2im{ Xi - x f )( Pl Q} + p f Qi 



(Here we have slightly changed notation, dropping the "0" suffices on the position and 
momentum parameters.) It can be verified that the exponential is a pure phase if and only 
Pf = Pi, Xf = Xi + Pit/m, in which case it becomes simply exp(itp 2 /2mh) . 

The initial and final conditions on the complex trajectories (J147P - (11481) are familiar from 
the literature, see in particular [20]. The semiclassical approximation (11461) is presented 
somewhat differently from from formulae in the literature, but it would seem to be equiv- 
alent. Our derivation, while maybe not as careful as previous derivations, is a substantial 
simplification. 

7 Concluding remarks 

The main results of this paper are as follows: After a detailed presentation of the complex 
WKB and BOMCA methods we showed how complex WKB can be derived from a saddle 
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point approximation in the path integral formulation of quantum mechanics. The path 
integral approach to the method explains the need to incorporate the contributions from 
multiple trajectories; the original formulation however is much more useful for practical 
applications, avoiding the cumbersome multiple integrals that arise when computing higher 
order correction terms from the path integral. In terms of methodology, the novel aspect of 
our path integral derivation was incorporation of the initial wave function into the integrand 
prior to computing the saddle points and the relevant behavior near them: It is this that gives 
rise to complex trajectories. Complex and real trajectory methods in quantum mechanics 
are complementary, not contradictory — complex trajectories are needed to propagate wave 
packet-type states as considered in this paper, whereas real trajectories suffice for WKB-type 
states. 

We then moved on to the path integral description of BOMCA. This required a further 
methodological innovation, the use of a general quadratic approximation in asymptotic 
analysis, as opposed to the standard Taylor approximation at the minimum. Using this 
more general asymptotic method we showed how to obtain BOMCA from the path integral 
(thus justifying the need for multiple trajectories in BOMCA too). In fact, from the path 
integral point of view at this stage BOMCA seems to be just one of many possible methods, 
a matter that merits further investigation. The overall picture of the relationship between 
complex WKB and BOMCA became clear. Both methods give rise to the same lowest order 
approximation to the wave function, the "classical wave function" (jSJ). In complex WKB 
this approximation is refined by multiplying the wave function by suitable factors of the 
form 1 + 0(h), while keeping the same trajectories x and their variations U. In BOMCA, it 
is the formula dHJ that remains the same, while 0(h) corrections are made to the trajectories 
x and the matrices U; these are changed depending on the order of the approximation. 

In section 6 we showed how our method of inserting the wave function into the path 
integral prior to making a saddle point approximation could be used to derive the coherent 
state propagator, measuring the overlap between an evolved Gaussian wave packet and 
another Gaussian state. Our derivation is a substantial simplifcation over previous ones. 
In the case of the coherent state propagator there is no alternative derivation to the path 
integral; for computations of the wave function, however, we emphasize that the derivation 
of the equations of complex WKB and BOMCA presented in section 2 is simpler than the 
path integral approach, which is only necessary to explain the need for multiple trajectories. 

There are a number of areas in which further work is necessary. First and foremost, this 
paper was intended to provide the theoretical backing for the numerical studies in [H [2j [3] , 
and having done this, we hope that further numerical studies will be undertaken, especially 
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in multiple dimensions. There are several areas in which more theoretical developments 
would be welcome. First, we have almost completely avoided in this paper any discussion 
of caustics (points at which the denominator in (jSJ) vanishes, rendering the approximation 
meaningless) and the related phenomenon of Stokes' lines. In the case of wave function 
approximations, the caustics and Stokes' lines are dependent on the time, and it is possible 
to write down equations describing their motion [59]. It is widely appreciated that the 
phenomena of caustics and Stokes' lines are "coordinate dependent", in the sense that 
they can be avoided (or moved) by working in momentum or phase space representations 
[531 EH [55l [56] . However, not enough has been done yet to make these ideas into efficient 
techniques for calculations. Strongly related to these questions is the more mathematical 
question of the nature of the singularities in the system of ODEs arising in complex WKB 
at a caustic, which we are currently investigating [57]. 

Another matter requiring further investigation is a better understanding of the (linear) 
decomposition of the wave function implied in (jSJ). Given an initial wave function is it 
possible (either abstractly or operationally) to write it as a sum of terms each of which 
evolves into one of the terms in the sum (jSJ)? Is the evolution by the Schrodinger equation 
or some other equation? Many ideas in these directions have been discussed by Poirier and 
collaborators [H M> 113 SSI M> HE] • 

Finally, we mention that we find the perturbed Newton's equations appearing in BOMCA 
to be fascinating. As mentioned above, from the path integral approach it emerges that the 
BOMCA equations are not unique, and we would like a way to select the version that 
emerges directly from the Schrodinger equation in section 2. We strongly suspect this to 
be related to some symmetry structure (recall that the underlying Newton's equations are 
Hamiltonian) , but have not yet found this structure. Understanding this might give us clues 
as to how to find nonperturbative BOMCA trajectories, that is trajectories that when used 
in (jSJ) give exact answers. In addition to these challenging, long-term goals, there is much 
to be done in seeking solutions of the first order BOMCA equations for specific systems 
and understanding, for example, the difference between the behavior of complex WKB and 
BOMCA near caustics. 
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Appendix A: Derivation of the classical wave function 
in the multidimensional case 

In this appendix we derive the multidimensional form of the classical wave function as 
described in the introduction and in section 3, see equations fl84|) . fl85|) . flHBT) and the text 
around them. The result differs slightly from standard results, specifically in the initial 
conditions obeyed by the classical paths (1841) and the function U flHBI) . and in any case the 
relevant calculations in the multidimensional case do not seem to have made it into most 
of the existing texts on path integration techniques, so we see fit to give at least the key 
details of the derivation. 



We start from the path integral in the form 

V(X,T) = JDxexp (I(£[x] + S init (x(0)))) 



where the integration is over all paths with x(T) = X. We use, in this appendix, an action 
of the form 






with a diagonal mass matrix; the case of a general mass matrix can be treated similarly. 
In the main text we quote results assuming all the masses rrii to be equal. We start by 
replacing x by x + e in the exponent and expanding to second order. Requiring the linear 
terms in e to vanish gives the classical equation of motion as well as the inital condition for 
x (1541) . The remaining terms give the approximation 

V>(X,T) « 5> i5M/ Vo(x(0) J De e Q[e] 

where 

Q ^ = ^{[ T ^ m Mtf-i2J2H(V) ij (^t))e i (t)e j (t)dt 

+ EE^(^ imt )^(x(o))5,(o)^(o) 
i=i j=i 

Here H(V) and H(S imt ) denote the matrices of second derivatives of V and S imt respectively. 
We proceed by discretizing the integrals in Q[e]. Bearing in mind that e(T) = and using 
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the trapezium rule we have 



fr(V r )« i (x(i))e i (t)e J -(*)dt 



^ff(V)y(x(0))e i (0)e i (0) 

T iV ~ 1 / /rT 

-^X>M«(*(^ 



rT\ /rT 
£i I — )e 



N J " J V N 

Using a forward difference approximation for the derivative of e(t) and a "leftbox" type 
approximation for the relevant integral gives the approximation 

2 

-2VJ 



T , x 2 , N ( , . 2 ^ /rT 



r=0 



A^ 



(Although this would appear to be a first order approximation, since both the methods for 
constructing the derivative and computing the integral are first order, it is actually second 
order; the first order errors in the methods exactly cancel each other.) Putting this all 
together gives: 

I M + F -M . . . \ 

-M 2M + F 1 -M 

-M 2M + F 2 -M 

-M 2M + F 3 -M 



Q[e] 



iN 
2~hT' 



:A 



y -M 2M + F N _! J 

Here A = (e(0) e(T/N) e(2T/N) . . .). Each entry in the matrix in the previous equa- 
tion is a d x d block matrix. M denotes the diagonal matrix with entries m,i, and the 
matrices F r are defined by 



( F o)ij 

\"r)ij 



^(5 init ),,(x(0)) - ^tf(V%(x(0)) , 
T 2 



-TriHWijlx 



rT\ 



N-l . 



N2 -y- nj y-y N J 

As a final step in the simplification of Q[e], we factor out factors of VM on the right and 
the left from each block in the above matrix. In this way we obtain 

I I + G -I \ 

-I 2I + G 1 -I 

-/ 2I + G 2 -I 

-I 2/ + G 3 -/ 



Q[e] 



iN 
2~hT 



AVM 



-I 21 + G 



'MA 



T 



N-l J 



(150) 
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where the matrix M. is a diagonal matrix with N copies of M on its main diagonal, and 
the matrices G r are defined by 



(G )ij 
[Lj r )ij 



^JWurfu \ N 



r F(5 ini %(x(0)) - ^#(V%-(x(0)) N 



1 T 2 



2N 2 ' 



^m-N^^V^KN 



N-l 



It remains to compute the determinant of the matrix in (11501) . To do this, we first apply 
block Gaussian elimination [62] to eliminate the blocks under the leading block diagonal, a 
process that does not affect the determinant. This gives a matrix of the form 



( P -I 

Pi -I 

P 2 -I 

P 3 -I 



\ 



P 



N-l J 



where 



Po — I + Co , 

Px = 27 + Gx-Po- 1 , 

P 2 = 2J + G 2 -Pr 1 , 

P/V-l = 2/ + (jjv-1 — P/v-2 • 

We wish to find det (P0P1P2 . . . P/v-i)- Define the matrices P r , r = 0, . . . , N — 1, by R r 
P PiP 2 . . . P r . Then we have 

Po = I + Gq , 

Px = (J + Go)(2/ + G 1 )-/ = / + 2Go + G 1 + G G 1 

and for 2 < r < iV - 1: 

P r = P r _2Pr— lPr = Pr— 2 (P — l(2i "I" Cj r J — i J = I\ r —iy21 -f- G>,) — P r _2 

Rewriting the last equation, for 2 < r < N — lwe have 

P r — 2P r _i + P r - 2 _ R ^_r* 

(7yJV)2 - ^-1 T 2 ^ ' 
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Taking the limit now as iV — > oo, we see that the sequence of matrices Rq, R\, . . . , Rn-i is 
replaced by a function R(t), defined by 

R(t) = R{t)G{t) , where G«(t) = — fffVWxft)) , 

supplemented with the initial conditions 

R(p) = I , %(0) = —^H(S™%(x(0)) . 

^mTifrTJ 

Finally, we write U(t) = R(t) T . Taking the transpose of all the equations above we see that 
U(t) satisfies 

Uit) = G(t)U(t) , U(0) = I, 4(0) = — ^i?(S' init ) ii (x(0)) . 

The determinant of the matrix in (11501) is simply det U(T), and (after checking the case of 
the free particle to fix the normalization) we deduce that 

^ Jdet U(T) 
as required. 
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